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Abstract 

Using random walk simulations we explore diffusive transport through monodisperse sphere 
packings over a range of packing fractions, cj), in the vicinity of the jamming transition at (f>c- 
Various diffusion properties are computed over several orders of magnitude in both time and 
packing pressure. Two well-separated regimes of normal, “Fickian” diffusion, where the mean 
squared displacement is linear in time, are observed. The first corresponds to diffusion inside 
individual spheres, while the latter is the long-time bulk diffusion. The intermediate anomalous 
diffusion regime and the long-time value of the diffusion coefficient are both shown to be controlled 
by particle contacts, which in turn depend on proximity to be- The time required to recover normal 
diffusion t* scales as (b — and the long-time diffusivity Doo ~ (0 — or Doo ~ 1/t*. 

It is shown that the distribution of mean first passage times associated with the escape of random 
walkers between neighboring particles controls both t* and Doo in the limit 4> ^ (fc- 
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Diffusive transport through heterogeneous media has been a long-standing area of in¬ 
terest [l, 2|, with wide-ranging applications in biology, geophysics and materials science. 
Heterogeneous media consisting of packed particles are of particular interest both as model 
systems and for industrial applications that involve powder processing. Analogous pro¬ 
cesses are electrical conduction in particle packs, relevant for energy storage applications 
(i.e. battery electrodes) e.g., j^; or molecular diffusion in porous materials with narrow 
inter-pore connections, relevant for biological applications (i.e. intercellular transport) jd]. 
Furthermore, the relationship between inhomogeneous material structure and bulk trans¬ 
port properties remains poorly understood in the context of multi-scale effects. Locally, 
material crystallline structure or interfacial effects can dominate transport, while at larger 
length and time scales, topology and connectivity prevail. Experimental investigations that 
simultaneously probe these length scales are challenging. In this work, we use detailed sim¬ 
ulations to explore the relationship between the bulk particle packing fraction 0, which can 
be measured and controlled experimentally; microstructural features that govern transport, 
such as the interparticle contact area; and diffusive transport over a broad range of length 
and time scales. 

A rich literature exists for effective medium theories of transport properties in hetero¬ 


geneous media 


lattices {d] as well as anomalous diffusion in disordered media 10|, lUl] from a stochastic 


, I 2 I, IsHSj , while much is known about diffusion on regular and disordered 


process point of view. In addition, work on percolation theory and in particn 
percolation has been directed at transport properties of random structures 


ar continuum 


12 


141. More 


recently, simulations of ballistic and diffusive 
overlapping spheres have been performed 


xansport using the Lorentz model of random 
16| . The results of that work are analyzed 


elegantly in a traditional continuum percolation-type framework (the so-called swiss cheese 
problem). Here, we are interested in conduction through the particle matrix and more im¬ 
portantly in systems that deviate from standard percolation (i.e. fractal) structures but 
that still show interesting scalings in the long-time limit (i.e. hyperuniform structures) 17 |. 

To this end, we focus on compressed, disordered packings of monodisperese spheres over 
a range of packing fractions 0, in the vicinity of the jamming transition 18| that takes 
place at a critical packing fraction 0c- 0 represents a convenient macroscale parameter that 
controls the nature of interparticle contacts: below 0c, packings lose mechanical stability, 
whereas for 0 > 0c, a number of scaling relations appear that determine various structural 
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FIG. 1. Disordered, jammed packings of monodisperse spheres: far from (pc (left; packing pressure 
= 4.25 X 10“^) and near pc (right; packing pressure = 4.29 x 10“^). Lines represent contacts 
between particles (not shown) with the thickness/shading indicating contact area (relative to the 
largest contact for that system). 


and mechanical properties of the packing jl9l-l2l|. Numerous studies j22| have highlighted 
the strong heterogeneity in the force or contact network. Such aspects of jammed packings 


lead to a number of unusual, or anomalous, mechanical and transport properties 23|, 


24|. 


Here, we take an explicit approach to exploring transport through jammed packings. We 
employ a stochastic technique based on a random walk algorithm 25| to track random walk¬ 
ers as they diffuse through the particle phase (equivalent to assuming a perfectly insulating 
surrounding medium), so that conduction cannot take place in the interstitial space. The 
interparticle contacts therefore play a crucial role in controlling the effective diffusivity, par¬ 
ticularly as 0 —)■ Pc, where the interparticle contact areas are small compared to the particle 
size. This key feature is illustrated in Fig. [H which shows two packings prepared at differ¬ 
ent packing pressures, or equivalently Ap = p — pc, where lines of thickness corresponding 
to the size of the contact areas join the centers of particles. For compressed systems far 
from the jamming threshold, not only is the average coordination number larger, but the 
contacts themselves have larger overlap areas. Closer to the jamming transition, contact 
areas decrease, and one would therefore expect diffusion through the packings to depend on 
proximity to the jamming transition. 

Indeed, we observe two distinct diffusion domains, one governed by diffusion of random 
walkers inside the particles, and the other by rare crossings of random walkers between 
neighboring particles. As p ^ pc, the separation in time scales for these domains diverges, 
and inter-particle hopping becomes the dominant, transport-limiting process. We therefore 
analyze the transport problem in the context of narrow escape/mean first passage time 26 |. 
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A recent analysis of diffusion from a sphere with multiple, well-separated, small absorbing 


windows 


27| provides the essential relationships between hrst passage times and interpar¬ 


ticle contact areas. Small changes in packing fraction, which affect interfacial contact areas 
between particles, ultimately lead to changes of several orders of magnitude in bulk prop¬ 
erties. Variations in interfacial areas in the current case of spheres with perfect contact 
are analogous to variations in interfacial contact resistance (in such systems, these could be 
controlled by system pressure, sintering processes, or particle surface treatments). Char¬ 
acterizing this contact resistance and connecting it to controllable bulk properties in the 
way we connect interfacial area to packing pressure is highly system-dependent; however, 
despite the idealized nature of the present case, these phenomena apply to a broad range of 
important systems. 

Our systems are composed of overcompressed, mechanically stable particle packings pre¬ 
pared using an established (de)compression protocol [28|. Each packing consists of N 
monodisperse, inelastic, frictionless spheres of diameter d and mass m in a cubic simula¬ 
tion box with periodic boundary conditions. All results presented here are in dimensionless 
units in which m = 1 and d = 1. Different packings were generated with particles interacting 
only on contact through a Hookean force law (spring force constant k = 1). By adjusting the 
packing fraction cj) over small intervals, multiple packings were generated at precise values 
of the packing pressure P, over the range 4 x 10' -8 < p 

< 4 X 10 For the systems studied 
here, we found 0c ~ 0.6383, which is consistent with the value of random close packing for 


frictionless spheres 


29[ |. Our major results are for N = 1000, but we also checked that system 


size did not affect our results using N = 10000, as shown in the Supplementary Material. 
One consequence of this packing generation protocol is the presence of rattlers - particles 
(or clusters of contacting particles) that have no contact neighbors. For typical jammed 
conhgurations close to 0c, less than 5% of particles are rattlers, which were removed from 
the packing prior to random walk simulations. 

We employed a random walker algorithm similar to that developed by Kim and Torquato 


25 


30 , 


3l| . The basis of the algorithm is a continous-time, off-lattice random walk within a 


heterogeneous structure, illustrated in the inset to Fig. [2l A key difference in our algorithm 
is that at each integration step, the time elapsed must be the same for all walkers. This 
leads to a slightly more complex simulation algorithm, but greatly simplihes the subsequent 
analysis of time-dependent properties. Specihcally, all random walkers must have the same 
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net displacement Axo (and therefore same elapsed time) at any given integration step [32 1. 
At the start of each simulation, ten times as many random walkers as spheres are randomly 
placed in the interior of particles. For a given walker at the start of a step, the distance to the 
edge of the particle in which it is currently located, denoted Axp, is first computed. If this 
distance is greater than Axo, the walker is translated a distance of Axo in a randomly chosen 
direction, and the time step is complete. The time associated with this move is Ax‘l/6Do, 
where Do is the bulk diffusion coefficient inside the particles, which we set to unity with no 
loss of generality. However, if Axi < Axp < Axo, where Axi is our minimal move resolution, 
the random walker is moved a distance Axp in a randomly chosen direction, but the time 
step is not completed. The process is repeated until the walker’s net displacement during 
the time step is exactly equal to Axo- This then completes the time step. If at some point 
Axp < Axi, the random walker is moved a distance of Axi in a random direction; if this 
move results in the random walker being located outside of any particle, it is reflected back 
into the interior of the particle by a distance Axi in a new random direction. The value 
of Axi can affect the measured properties, and must be varied accordingly to resolve the 
smallest relevant length scale in the problem, which here is the size of the contact area (i.e. 
smaller A0 require smaller Ax*); tests of convergence with respect to Axi are shown in the 
Supplementary Material. 

Random walk simulations were carried out on the different packings for a wide range of 
pressures. Figures [2] and [3] show key measures characterizing the diffusion of random walkers 
in jammed particle packs for several values of A0. 





FIG. 2. Mean squared displacement (MSD) of random walkers in jammed packings spanning 
10-6 < Acj) <10 ^. Insets: The anomalous diffusion exponent a{t) and a schematic of the random 
walker simulation. Color legend is the same as in Figure [3l 
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FIG. 3. Evolution of the effective diffusion coefficient for the different packings. 

As seen from the plot of mean squared displacement (MSD) in Fig. [21 the short-time 
diffusion behavior is Fickian (MSD ~ t) for all A0. This corresponds to diffusion of random 
walkers inside individual particles, on time scales shorter than that at which random walkers 
encounter the particle boundary. At intermediate times (starting at t 10 ^), a subdiffu- 
sive plateau develops, which becomes broader as 0 ^ (l)c- At sufficiently long times, linear 
diffusive behavior is once again recovered for all systems as the walkers explore the entire 
packing. The time required to reach the long-time linear regime increases with decreasing 
A0. This time-dependent diffusion behavior can be expressed either as a time-dependent ef¬ 
fective diffusion coefficient D(f), such that MSD = D{t)t/6, or a time-dependent anomalous 
diffusion exponent a{t), MSD ~ /Q. The inset of Fig. [2] shows a{t) = dlogMSD/dlogf 

for several values of A0. The Fickian diffusive regimes at short and long times correspond 
to a ~ 1. Alternatively, Fig. |3] shows the time-dependent behavior of D{t). At short times, 
D{t) Do, corresponding to the diffusion coefficient inside the particles. At longer times, 
D{t) decreases due to the trapping of random walkers inside particles, which becomes more 
pronounced as 0 ^ 0c due to smaller interparticle contacts. The long-time plateau value 
corresponds to diffusion at time scales beyond those associated with interparticle hopping, 
and thus reflects the bulk homogeneous limit. 

In order to quantify the transition between the various diffusion regimes, we dehne the 
characteristic time such that «(to 95 ) = 0-95 near the transition from the intermediate 
subdiffusive regime to the Fickian regime. This represents the characteristic time at which 
bulk Fickian diffusion is recovered, tg.gs is plotted in Figure |4|(a) as a function of A0. A 
clear scaling of t* ~ (0 — 0c) is observed, with perhaps the exception of large A0 values. 
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FIG. 4. (a) Characteristic time gg at which Fickian diffusion is recovered. Dashed black line 
corresponds to (b) Late-time value of the diffusion coefficient D^o as a 

function of Acp. Line corresponds to Dqo = 2.0(0 — 


which lie outside of the traditional jamming-scaling regime 20|]. Similarly, in Fig|4] (b), the 
long-time plateau value of the diffusion coefficient (late-time values in Fig [3]) is shown as a 
function of A0. Note a clear scaling of Doo ~ (0 — 0c) 

To understand the scaling behavior observed for t* and Doo, we analyze the diffusion of 
random walkers in the context of a narrow escape, or first passage time 2^. Here the first 
passage time is the time for a random walker to escape the interior of a given particle to 
any of its neighboring particles. To escape, a random walker must reach the small circular 
contact areas that separate a particle from its neighbors. The analogous problem has been 
studied by Cheviakov et al |27| , who derived an asymptotic analytical solution for the mean 
first passage time MFPT(x) as a function of the starting position x and the size and relative 
orientation of small circular openings on the surface of a spherical particle. The analytical 
solution requires that the radius of the circular opening's be much smaller than the radius 
of the particle r and that they are well separated 27|. These conditions are satisfied for 
all cases of interest here. For a given particle, we compute the volume-averaged MFPT, i.e. 
the MFPT averaged over all starting positions x and random walk realizations inside the 
particle (for a full description, see the work of Cheviakov et al j^, in particular eq. 2.44 of 
that work). We refer to the volume-averaged MFPT as i. 

Based on the structures of the packing configurations, we have computed t for A0 = 
9.8 X 10“^, and have compared this result to our random walk data in Fig. E](a). The 
agreement is very good; small discrepanices are due to convergence errors in the simulation 
data (for these simulations, convergence errors are much larger than other data previously 
discussed). Given this satisfactory correspondence, we use the analytical result of Cheviakov 


et al 


27j in our subsequent analysis. For each system, the median value of t, denoted tm, is 
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FIG. 5. (a) Comparison of volume-averaged MFPTs computed from simulations (red) and from 
the analytical expression 0 (black) for A(/> = 9.8 x 10 The main fignre shows the comparison 
of the histograms, the inset shows the direct comparison of values for all particles. The solid line 
in the inset has a slope of unity and zero intercept, (b) Relationship between tg^/tm and MFPT 
for all A(j). Inset shows corresponding MSD at t = tg^ 


designated as the system average MFPT. In Figure [5](b), we plot as a function of A0. 

At all but the highest A0 values, the ratio is constant, indicating that the characteristic 
time t*, which signals the crossover from the subdiffusive to the long-time, bulk diffusive 
behavior is set by the MFPT, i.e. t* ~ Hence, large-scale diffusion, as 0 ^ is limited 
by the ability of a random walkers to escape from an average particle to its neighbor, which 
in turn is governed by the average size of the interparticle contacts. 


The constant value of ~ 10 to which the tgs/tm ratio converges in Figure [5](b) is not 
readily explicated, since it is convoluted by the fact that is a median value of the volume- 
averaged MFPT. Instead, we plot the MSD of the random walkers at t = tg^ in the inset 
of Fig|Sl(b). This quantity, which we denote MSD*, converges to a constant value of ~6 as 
A0 decreases. This can be readily interpreted as the average distance, \/MSD* ~ 2.5, that 
a random walker must travel before its trajectory exhibits regular diffusion characteristics. 
In other words, in the limit 0 (p^ where diffusion is dominated by rare event interparticle 

hopping, ~2-3 such hops are needed to reach the diffusive regime. Thus, the anomalous 
diffusive behavior is localised on the particle scale and not the direct result of any large 
length scale phenomenon. It is interesting to note, however, that this is consistent with the 
distance at which the radial distribution function, g{r), loses history of the local packing 
structure and is likely consistent with the idea that the directionality of the discrete contact 
angles between particles is averaged out. 


We have established that t* ~ 
leading order, the analytical result 


tm, and Fig. |4] shows that t 


27( 1 yields t ~ and therefore tr. 


(0-0,)-o-^ To 
where Tc is 

















the radius of the contact region. By geometric arguments (0 — 0c) ~ <5 for small 6, where 6 
is the overlap between contacting particles i and j, 6 = d — Vij, and rij is their center-center 
separation. Similarly, Vc ~ 5^'^ for small 6. Combining all of these, we hnd: 


t* r. 


(0-0C 


\-0.5 


( 1 ) 


which is consistent with the results of Fig. ID^a) 33|. 

The scaling behavior of the long-time effective diffusion coefficient Dao can also be readily 
explained by a similar argument. In the long-time limit, normal diffusion is recovered, 
such that MSD ~ t. For any t* value selected so that a is arbitrarily close to unity, the 
corresponding MSD* will be constant at low A0 values (see inset of Figure [5](b) for tg.gs)- 
The long-time diffusion coefficient in this case is then simply Doo = MSD*/6f*, and so 
Doo ~ 1/t*. Since t* ^ {(j) — 0c) it follows that D^o ~ (0 — 0c)the scaling behavior 
noted in Figure |l](b). 

In summary, we have presented random walk simulations of diffusive transport in dis¬ 
ordered particle packings in the vicinity of the jamming transition. Two distinct linear 
diffusive regimes where MSD ~ t are observed at short and long time scales, separated by 
an intermediate subdiffusive regime. The first regime corresponds to regular diffusion in 
the interior of particles, while the second regime is the expected long-time bulk diffusion 
behavior. The time scale separation between the two regimes grows as the volume fraction 
approaches the critical jamming value. This is a result of the corresponding decrease in the 
size of interparticle contacts, which leads to the average diffusion behavior being controlled 
by the ability of random walkers to escape from one particle to neighboring coordinated 
particles, so-called narrow escape. This behavior can be analyzed in the context of a mean 
hrst passage time. We observe scaling of t* (0 — 0c) and Doo ~ (0 — 0c)°'^, where t* is 
a characteristic time at which linear diffusion is recovered and Doo is the long-time effective 
bulk diffusion coefficient. Since nearest-neighbor inter-particle contacts control the large- 
scale transport properties in the particle packs considered here, this behavior is not unique 
to systems near the jamming point since it is not the presence of disorder that controls the 
anomalous diffusive behavior but rather the size of the particle contacts themselves. There¬ 
fore, we expect [3^ the qualitative trends, as well as the scaling relations reported here, to 
be independent of the interparticle potential as well as packing geometry. It is hoped that 
this work will motivate detailed experimental studies of tranport dynamics at intermediate 
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scales in systems dominated by interfacial resistance, e.g. pulse thermography 
particle packs. 
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